/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		triangle.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date:  2020-01-17 16:52
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-08-04 01:26
*
*   Description:
*
=============================================================================*/
#ifndef I_TRIANGLE_H
#define I_TRIANGLE_H

#include "polygon_interface.hpp"
#include "triangle_shape_function.hpp"

namespace geoxfem {
namespace geometry {

template <typename DimensionType, typename degree>
struct Triangle : PolygonInterface<DimensionType, degree> {
  DEFINE_KEY(Triangle)

  static constexpr size_t dimension = traits::dimension<DimensionType>::value;

  using Shape      = ShapeInterface<DimensionType, degree>;
  using Polygon    = PolygonInterface<DimensionType, degree>;
  using point_type = typename traits::point_type<DimensionType>::type;
  using Gauss      = model::GaussPoint<DimensionType>;

  /*constructors and destructor{{{*/

  /**
   * @brief constructor.
   *
   * @param _listPoint
   * @param _pMesh
   */
  Triangle(const std::vector<point_type>& _listPoint) : Polygon{_listPoint} {}
  Triangle()          = default;
  virtual ~Triangle() = default;
  /*}}}*/
  /*misc{{{*/

  /**
   * @brief getNumberOfPoint
   *
   * @return
   */
  size_t getNumberOfPoint() const final override {
    return degree::value * 3;
  }
  /**
   * @brief get Class Name
   *
   * @return
   */
  std::string getClassName() const final override {
    return std::string("Triangle") + Polygon::poly_degree[degree::value]
           + std::to_string(dimension) + "D";
  }
  /*}}}*/
  /*shape function using reference coordinates{{{*/

  /** calculate standard shape function using local coordinates
   *
   * @param lcoord
   * @return
   */
  VectorXR getStdNByLocal(const VectorXR& lcoord) const final override {
    return triangular_shape_function<dimension, degree::value>::getStdNByLocal(
      lcoord, *this);
  }
  /** calculate standard shape function derivation dNdlx using local
   * coordinates
   *
   * @param lcoord
   * @return
   */
  MatrixXR getStddNdlxByLocal(const VectorXR& lcoord) const final override {
    return triangular_shape_function<dimension,
                                     degree::value>::getStddNdlxByLocal(lcoord,
                                                                        *this);
  }
  /*}}}*/
  /*shape function using global coordinates{{{*/

  /** also you should get shape function N with global coordinate gcoord
   *
   * @param gcoord
   * @return
   */
  VectorXR getStdNByGlobal(const VectorXR& gcoord) const final override {
    return triangular_shape_function<dimension, degree::value>::getStdNByGlobal(
      gcoord, *this);
  }
  /** then global coordinate can be yield by function above
   *
   * @param gcoord
   * @return
   */
  VectorXR global2local(const VectorXR& gcoord) const final override {
    return triangular_shape_function<dimension, degree::value>::global2local(
      gcoord, *this);
  }
  /*}}}*/
  /*gauss point related{{{*/

  void getStdGaussPoints(std::vector<Gauss>& gp_list) final override {
    size_t   num_gp = model::i_default_gauss_num::set_get().at("triangle");
    MatrixXR lcoord(dimension, num_gp);
    VectorXR lweight(num_gp);

    switch(num_gp) {
    case 1:
      lcoord << 1. / 3., 1. / 3.;
      lweight << .5;
      break;
    case 3:
      lcoord << .5, .5, 0, .5, 0, .5;
      lweight << 1. / 6., 1. / 6., 1. / 6.;
      break;
    case 6:
      lcoord << 0.81685, 0.09158, 0.09158, 0.10810, 0.44595, 0.44595,  //
        0.09158, 0.81685, 0.09158, 0.44595, 0.10810, 0.44595;
      lweight << 0.05498, 0.05498, 0.05498, 0.11169, 0.11169, 0.11169;
      break;
    case 7:
      lcoord << 0.33333, 0.79743, 0.10129, 0.10129, 0.47014, 0.05972, 0.47014,
        0.33333, 0.10129, 0.79743, 0.10129, 0.05972, 0.47014, 0.47014;
      lweight << 0.11250, 0.06297, 0.06297, 0.06297, 0.06620, 0.06620, 0.06620;
      break;
    case 12:
      lcoord << 0.87382, 0.06309, 0.06309, 0.50143, 0.24929, 0.24929, 0.63650,
        0.63650, 0.31035, 0.31035, 0.05315, 0.05315,  //
        0.06309, 0.87382, 0.06309, 0.24929, 0.50143, 0.24929, 0.31035, 0.05315,
        0.63650, 0.05315, 0.63650, 0.31035;
      lweight << 0.02542, 0.02542, 0.02542, 0.05839, 0.05839, 0.05839, 0.04143,
        0.04143, 0.04143, 0.04143, 0.04143, 0.04143;
      break;
    case 16:
      lcoord << 0.33333, 0.65886, 0.17057, 0.17057, 0.89891, 0.05055, 0.05055,
        0.08141, 0.45929, 0.45929, 0.00839, 0.00839, 0.26311, 0.26311, 0.72849,
        0.72849,  //
        0.33333, 0.17057, 0.65886, 0.17057, 0.05055, 0.89891, 0.05055, 0.45929,
        0.08141, 0.45929, 0.26311, 0.72849, 0.00839, 0.72849, 0.00839, 0.26311;
      lweight << 0.07216, 0.05161, 0.05161, 0.05161, 0.01623, 0.01623, 0.01623,
        0.04755, 0.04755, 0.04755, 0.01362, 0.01362, 0.01362, 0.01362, 0.01362,
        0.01362;
      break;
    default:
      throw "illegal gauss point number!";
      break;
    }
    for(size_t i = 0; i != num_gp; ++i) {
      VectorXR gcoord = this->local2global(lcoord.col(i));
      _real    gweight =
        lweight[i] * this->getJacobByLocal(lcoord.col(i)).determinant();
      gp_list.emplace_back(Gauss{gcoord, gweight});
    }
    return;
  }

  /*}}}*/
  /*io related{{{*/
  /**
   * @brief getVTKTypeID
   *
   * @return
   */
  size_t getVTKTypeID() const final override {
    return getVTKTypeIDImpl<degree>();
  }
  template <typename d>
  size_t getVTKTypeIDImpl() const {
    return 53;
  }
  template <>
  size_t getVTKTypeIDImpl<mpl1>() const {
    return 5;
  }
  template <>
  size_t getVTKTypeIDImpl<mpl2>() const {
    return 22;
  }

  std::string getVTKPointOrder() const final override {
    return get_vtk_point_order_impl<degree>();
  }
  template <typename d>
  std::string get_vtk_point_order_impl() const {
    throw "not implemented yet";
  }
  template <>
  std::string get_vtk_point_order_impl<mpl1>() const {
    return Shape::getVTKPointOrder();
  }
  template <>
  std::string get_vtk_point_order_impl<mpl2>() const {
    /* Create an output string stream */
    std::ostringstream streamObj;
    /* Add double to stream */
    streamObj << this->getPoint(0).getIndex() << " "  //
              << this->getPoint(4).getIndex() << " "  //
              << this->getPoint(2).getIndex() << " "  //
              << this->getPoint(5).getIndex() << " "  //
              << this->getPoint(3).getIndex() << " "  //
              << this->getPoint(1).getIndex();
    /* Get string from output string stream */
    return streamObj.str();
  }

  /*}}}*/
};

}  // namespace geometry
}  // namespace geoxfem

GEOXFEM_REGISTER_TEMPLATIZED(geoxfem::geometry::Triangle);

#endif /* I_TRIANGLE_H */
